Emergence of Critical Phenomena in Full Configuration Interaction Quantum Monte 

Carlo 
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There has been recent literature discussion on the origin and severity of the 'sign problem' in 
full configuration interaction quantum Monte Carlo (FCIQMC) and its 'initiator' adaptation (i- 
FCIQMC), methods of interest and potential because they allow for exact (FCI) ground-state so- 
lutions to be obtained often at a much reduced computational cost. In this study we aim to use a 
simple order parameter, describing the 'sign structure' of the stochastic wavefunction representation, 
to empirically characterise the fundamentally different collective behaviour of the walker population 
in both methods. 



Full configuration interaction quantum Monte Carlo 
(FCIQMC) PJ, and its initiator adaptation[l Q (i- 
FCIQMC), are electronic structure methods which at- 
tempt to itcratively solve the imaginary-time Schrodingcr 
equation within a stochastic dynamic of walkers spanning 
a Slater determinant space. This approach has a great 
deal of potential because it has been shown to effectively 
diagonalize extremely large spaces and obtain ground- 
state ener gies, in molecules and a number of model 
systems (a 1 1 0Ml~j | . 

Although this method is proving effective in a range 
of problems, analysis of why it is so successful at re- 
ducing the computational cost relative to FCI is still 
a topic of discussion 0, HI, [l3|-[l6|. In this communica- 
tion, we investigate the abrupt emergence of the 'sign 
structure' of the wavefunction in FCIQMC to the exact 
solution by making an analogy with symmetry-breaking 
phase transitions and related critical point phenomena. 
This is contrasted with the smoother onset of order in 
the i-FCIQMC method. In the spirit of previous work 
analysing the so-called 'sign problems' of other quantum 
Monte Carlo methods, it is hoped that work such as this 
to understand its origin and amelioration will lead to 
better strategies and lower cost in its treatment (l7l - [30| . 

The crux of the FCIQMC method is that the finite- 
basis Schrodinger equation (H^o = E^o; ^>q = 
^;Ci|£>i)), is recast as a set of master equations gov- 
erning the dynamics of a walker population in Slater de- 
terminant space, 
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where the aim is to find ^p 1 = and recover S as the 
correlation energy. The walkers themselves take a weight 
of ±l[lH and in a converged calculation it is intended 
that they be distributed such that on average, Cj oc N\. 
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These dynamics consist of different stochastically real- 
ized events, but in particular positively and negatively 
signed walkers are removed from the simulation by an 
'annihilation' criterion. Walkers of opposite signs simul- 
taneously on the same site are removed at the end of 
each iteration, this procedure ensuring that each deter- 
minant has a definite sign once enough walkers have been 
introduced into the simulation. This greatly ameliorates 
the FCIQMC 'sign problem', which has been investigated 
by Spencer, Blunt and Foulkes[l3j, and is responsible for 
the required complexity of the instantaneous stochastic 
solution. They isolated a solution which the FCIQMC al- 
gorithm without sufficient annihilation is unstable with 
respect to, characterised by the ground state of a modi- 
fied Hamiltonian matrix, H'^ = — for i 7^ j. 

This explained the observation of a so-called annihi- 
lation plateau in FCIQMC: a 'stalling' in population 
growth at a characteristic walker number (N c ). Above 
this walker number, the ground-state energy for the true 
Hamiltonian can typically be extracted, and below this 
walker number, the solution is contaminated by the afore- 
mentioned state. 

In the initiator adaptation to FCIQMC, an additional 
criterion is added to the step where walkers are created 
(spawning). The space is divided into those exceeding a 
certainly system-dependent parameter n a dd, termed 'ini- 
tiators'. Spawning is conducted such that is zeroed 
when spawning onto a site with no walkers, if the origin 
of the spawning attempt was not an 'initiator'. The effect 
of this is to dynamically modify the Hamiltonian used to 
generate "J (r + St) from '5 (t) . This can introduce an 
error at lower walker number due to the truncation of 
the instantaneously available space. Much of this error 
will be removed in the time-averaging, however system- 
atic error can remain, and is referred to as the 'initiator 
error'. 

An FCIQMC or i-FCIQMC simulation in general be- 
gins with a single, positively-signed walker on a reference 
determinant (typically the HF determinant) [32j . From 
this point, the population is grown by setting the shift 
parameter to be fixed such that S > E . Typically, at 
some point in the simulation the population is relaxed 
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from this fixed-shift mode, and the shift is updated in a 
self-consistent fashion (33| . 

We chose as our model system the beryllium dimcr 
at a near-equilibrium geometry of 3.8329 a.u., in a cc- 
pVTZ basis set, with frozen core electrons. Our molec- 
ular integrals were calculated by the QCHEM package^ 
(-Ehf = —29.1123 a.u.). This was chosen predominantly 
because it is a system with a small determinant space 
of Npci = 346 485, with a comparatively large plateau 
height of A c ~ 164 000 walkers (there is a narrow gaus- 
sian distribution of A c between different simulations, 
S'=0.0 a.u., with a time step of 0.001 a.u.). It is also 
an example of a system trivially solved by j-FCIQMC 
with very little initiator error at modest walker numbers. 
A typical simulation profile for FCIQMC, showing the 
growth of the walker population, is shown in Fig. Hal At 
approximately r = 3, the population enters the so-called 
annihilation plateau and stabilises for a certain length of 
time before exiting (around r = 50). The length of time 
spent in the plateau varies between simulations and has 
an exponential distribution in its tail. No such plateau 
exists in an equivalent initiator simulation for Be2 and 
for all molecular systems studied to date. 

We seek a single parameter which we can use to moni- 
tor the convergence of the instantaneous stochastic wave- 
function representation, 



*(r) = £a|A> 
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to the true ground-state (FCI) wavefunction 

*o = E c iVi>- ( 3 ) 

or the degenerate alternative —^o, where c^ ' and C are 
normalised. Motivated by discussion of 'sign structure' 
in previous papers we construct the following parameter, 
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where N^t is the number of non-zero weighted determi- 
nants in and the prime above the sum indicates that 
the sum should not be taken over any component where 
the numerator is zero. This expresses the proportion of 
the signs that are correctly-aligned with the final wave- 
function, with the following interpretation: 

= 1, signs completely aligned with + 
> 0, signs more aligned with + than — 
( =0, randomly aligned signs 

< 0, signs more aligned with — than Wo 
= — 1, signs completely aligned with — Wo- 

(5) 

Due to the similarity with order parameters in many 
branches of statistical physics we can also consider this 



as such. We will also employ an occupied- space order pa- 
rameter in this study, where the parameter is re- weighted 
by Adot/A occ , where A occ is the number of non-zero 
weighted determinants in (r) . We will simply denote 
this jf^m. This is in general larger than m (although 
has the same limiting values of ±1) and has the effect of 
only finding the ratio of sign alignment with in the space 
which has walkers in it, rather than over the whole space, 
and in particular has a starting-point of ±1 rather than 
± jj— . We will also compare these with an overlap value 

c^ * 1 • C, which includes consideration of the magnitude 
as well as sign overlap of the wavefunctions. 

For FCIQMC simulations on the beryllium dimer di- 
rectly corresponding to Fig. Ila[ calculated order parame- 
ters are shown in Fig. Ilbl with respect to different walker 
populations, A w . We take the approach of only showing 
the instantaneous order parameter of a single trajectory 
at a time, rather than attempting to deal with stochas- 
tic error, noting that these are extremely similar to data 
from averages taken by going into so-called variable shift 
mode. This is because the equilibration time for this 
system is close to zero and the qualitative behaviour at 
each point insensitive to the starting random-number se- 
quence. This situation is discussed extensively in another 



paper by the authors [11 



The rising magnitude of the order parameter results 
from gradual rise in the overlap between the simulation 
wavefunction and exact wavefunction (c*- ) -C, also shown 
in Fig, lib"]) . In particular, the graph shows two regions. 

For A w < A c , the order parameters fluctuate about 
zero. This represents no average net alignment of the 
signs, or overlap, with either of +^0 or —Wo. In pre- 
plateau simulations, or simulations without sufficient an- 
nihilation, the ground-state solution is contaminated by 
a solution from an effective Hamiltonian (H~ = — \H{j\ 
as described earlier (HI)- Since the effective Hamiltonian 
only has negative off-diagonal matrix elements, the resul- 
tant eigenvector has same signed coefficients, although 
there is still a degenerate solution corresponding to a 
global sign change. The bulk of the positively (nega- 
tively) signed walkers are distributed according to the 
positively (negatively) signed solution in this phase, re- 
sulting in no net instantaneous overlap with ±Wo. This is 
the region in which the sign of the reference determinant 
can be seen to flip in some svstems[14|. 

For A w > A c , the order parameters become non-zero, 
showing significant emergence of either —Wo or +Wo as 
the dominant solution (we have shown an example of 
both types of trajectory). The growth in this parameter 
m is then apparently logarithmic in walker number. This 
contrasts the true overlap which grows extremely rapidly. 
The rapid growth in this overlap causes the projected 
energy and shift estimators for the correlation energy to 
be exact beyond exit from the plateau [H-[4[. 

The discontinuity at N w = N c corresponds directly to 
entry into the plateau region. In fact, the overlap order 
parameter c^ ) ■ C grows to ^0.6 before the a clear exit 
from the plateau is seen, a value which fluctuates strongly 



3 



10"' 



FCIQMC 
i-FCIQMC 



10" J 



10" 



10 1 



10^ 



Imaginary time elapsed r (a.u.) 



1.0 



0.5 



U ■ 
O 



0.0 



-0.5 



-1.0 



m 


r 










— c<°).C 




N W =N C 




; 

: 

: 

; 





10* 10 3 10" 10' 

Population size N w (io 5 walkers) 




Rescaled effective temperature : 



(a) (b) (c) 

FIG. 1: (a) In FCIQMC, there naturally arises a plateau in walker growth termed the annihilation plateau 
(iV c = N w ~ 164, 000). This feature is absent in i-FCIQMC. (b) Above this critical walker number, there is a change 

in the order parameter from oscillating about zero to having a non-zero value. The sign is trajectory-dependent 
(dashed lines correspond to the same parameters as the solid lines, with a different trajectory). The speed at which 
this reaches ±1 is dependent on the order parameter (discussed in the text), (c) Transforming 1/N W to an effective 

temperature (T), the critical exponent for small values of the order parameter are found to be P = 0.77 ± 0.03. 



with iV w and imaginary time. 

The slow approach of the order parameter m to unity 
is due to the instantaneous wavefunction bearing little re- 
semblance to the final wavefunction "Jo for low-amplitude 
determinants, only being on average able to yield these 
coefficients 0- Stochastic fluctuations therefore not only 
modify the coefficients of the eigenvector but also the sign 
structure. In principle, we could average the wavefunc- 
tion over imaginary time, whereupon we would expect 
the average to approach ho- 
using our physical order parameter analogy, we can 
think of the walker population N w as being related to 
an effective temperature scale T = 1/ (N w ). The plateau 
height N c can be seen to be a critical temperature T c , 
and at this critical temperature the population settles 
into an ordered phase of either +^0 or — ^o- Ergodicity 
is broken and the other ordered phase is not explored 
beyond entry into the plateau. This phase transition, 
characterised by our order parameter, can be seen to be 
second-order with a critical exponent defined by, 



with p ~= 0.77 ± 0.03 (Fig. [TcJ). The region of fitting 
is the linear region taken by the distribution of points 
beyond 10 -2 , and varying the plateau height by ±2 000 
(from N c = 164 000) causes deviation in our estimate of (3 
by at most ±0.05. This /3 is consistent with the data from 
overlap order parameter c' ) ■ C, and although cannot be 
independently found from a free-fit due to noise in the 
data demonstrates that measurement of P is consistent 
between different measures of this ordering. 

There remains much to be understood about the mean- 
ing of this p parameter and a complete investigation of 
this is well beyond the scope of this initial study. We 



speculate that its value being greater than a (Landau) 
mean-field value of P = | means that there is anti- 
cooperativity between the presently ordered phase and 
the surrounding space; a new determinant needs to be 
consistently ordered with respect to a growing set of or- 
dered determinants. We expect it would be possible to 
link this to behaviour in nature, and note that this value 
is similar to that observed by simulations of frustrated 
spin- glasses [35]. Although universal between different 
methods of measuring the ordering, we would expect it 
to be strongly system-dependent and hence might be a 
measure of how intrinsically entangled the wavefunction 
is. 

Having shown that our order parameter m can pro- 
duce the anticipated behaviour for FCIQMC, we turn 
our attention to the initiator adaptation of FCIQMC 
(i-FCIQMC). We begin by noting that Be2 is almost 
immediately solved for almost any walker number in i- 
FCIQMC. The dimer is actually an extreme example of 
a difficult system for FCIQMC, and a straight-forward 
system for z-FCIQMC. 

In Fig. I2a[ the same order parameter analysis is re- 
peated for the initiator approximation. Here, averaged 
values obtained in variable shift mode are shown by solid 
circle on the plot, and the instantaneous and averaged 
values overlay one another with a slight shift in N w . This 
mismatch is anticipated, and is due to S ^ E COTr as de- 
scribed in Ref. [ll|. 

The order parameters for i-FCIQMC are very different 
to those for FCIQMC. The whole-space order parame- 
ter m, rises continuously remaining positive throughout. 
Correspondingly, there is no sign-flipping of the refer- 
ence determinant as in pre-plateau FCIQMC simulations 
for this system. The occupied-space order parameter, 
jVdet m ^ gj so remams positive for the same reason. How- 
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FIG. 2: (a) In i-FCIQMC, there is no annihilation plateau and the order parameter does not change sign, remaining 
at the same sign as the original walker distribution (positive in this case). There is immediate monotonic growth of 
the whole-space order parameter m. In contrast, the reduced-space order parameter initially falls, due to increasing 

distribution of the walker population across the space, before rising as this is overtaken by these signs aligning. 
Filled circles represent averaged values taken by going into variable shift mode, whereas the dashed line traces out 
the instantaneous values during walker growth, (b) When the order parameter is plotted with respect to the 
population on the Hartree-Fock determinant, FCIQMC and i-FCIQMC are much more similar, (c) For a given 
population size, there is greater annihilation in i-FCIQMC than FCIQMC. 



ever, it falls from its initial starting point of +1 (or in gen- 
eral ±1), as the walkers distribute themselves across the 
space at a quicker rate than sign coherence can spread. 
The latter then catches up as more of the space becomes 
aligned. The two order parameters m and -j^^m differ 
significantly in i-FCIQMC whilst they do not differ by 
much in FCIQMC. This is in part because the occupied 
space in i-FCIQMC is typically (always for Be2) smaller 
than the occupied space for FCIQMC. 

Finally, Fig. I2bl shows that, from the point of view of 
the population at the Hartree-Fock determinant, there 
is a remarkable similarity between order parameters 
in FCIQMC and i-FCIQMC. This has been observed 
previously as the idea that the HF determinant only 
grows its population in the post-plateau phase of the 
calculation [36(, and further relates the symmetry-broken 
phase of FCIQMC to the whole of an i-FCIQMC simu- 
lation. 

These data are consistent with the idea that the annihi- 
lation plateau is removed by the 'initiator' adaptation for 
systems in which it is effective, which can also been seen 
in enhanced annihilation rates in i-FCIQMC compared 
with FCIQMC (Fig. [5c]). We agree with Spencer et aZ.0 
that the restriction of the simulation to a sub-space im- 
posed by linking VP (r) and H-^ (r) leads to increased anni- 
hilation because the walkers are able to 'meet' each other 
more. On the other hand, there is no non-trivial link be- 
tween this and a sub-space diagonalization, since energies 
are much-improved compared with those that are pro- 
duced by a diagonalization in the occupied subspaceQ- 
Looking at the beginning of the simulation, there is much 
slower growth in i-FCIQMC than in FCIQMC (Fig.©, 
which is because walkers are being admitted at a much 



slower rate into the simulation. This restriction, and 
the smooth growth in the order parameter, happens at a 
much lower walker number than significant annihilation 
begins. 

In summary, we have introduced a parameter (m) 
which measures the extent of agreement between the 
signs of the stochastic wavefunction representations in 
FCIQMC and i-FCIQMC and the true ground-state so- 
lution. When appropriately normalised, this measures 
±1 for perfect sign agreement with In FCIQMC 

simulations of the beryllium dimer, m fluctuates about 
zero for small populations, taking both positive and neg- 
ative values as the sign of the reference determinant (and 
all determinants) flips. After the annihilation plateau is 
entered, for a single trajectory the sign of m is 'locked 
in' to either approach +1 or —1 as the walker population 
is grown. This is linked to a phase transition at an ef- 
fective temperature given by an inverse walker number. 
In contrast, i-FCIQMC has no annihilation plateau, sign 
flipping, or observed phase transition behaviour in the 
order parameters, instead immediately growing towards 
+1. We relate to increased annihilation in the system 
due to the dynamical Hamiltonian, which contrains the 
instantaneous space and enhances the removal of the sign 
problem in i-FCIQMC by annihilation. 

As with any empirical or numerical study of such phe- 
nomena, one important issue to address is the trans- 
ferability of these findings. In particular, Be2 is a sys- 
tem with a pronounced plateau in FCIQMC, but that is 
largely problem- free in i-FCIQMC. We have effectively, 
therefore, used it to show how we believe i-FCIQMC to 
be working when it works at its best. In particular, these 
conclusions would probably not be true in an i-FCIQMC 
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